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Generalized  Multistep  Methods  and  Applications  to 
Satellite  Orbit  Trajectory  Computation 

by 

James  Dyer 
ABSTRACT 

Recent  ideas  in  the  theory  of  multistep  methods  of  solving  a  differential 
equation  of  the  first  order  are  extended  to  methods  of  solving  the  special 
second  order  equation  y"  «  f(x,y).  A  slight  modification  of  the  usual  multi- 
step  method  permits  £  significantly  higher  order  approximation  of  the  difference 
equation  to  the  differential  equation  without  loss  of  stability. 

The  method  of  constructing  the  generalized  difference  equations  is  based 
on  a  quasi-Hermite  polynomial  approximation.  An  outline  of  this  theory  is 
given  along  with  some  related  unsolved  problems.  This  method  permits  the  con¬ 
struction  of  new  classes  of  stable  difference  equations  with  high  order  of 
accuracy  for  solving  both  a  first  order  differential  equation  and  the  above 
special  second  order  equation. 

Some  of  the  new  methods  have  been  tested  in  experiments  including  the 
computation  of  an  unperturbed  satellite  orbit  trajectory.  Machine  time  used 
and  accuracy  obtained  are  compared  with  a  standard  multistep  method. 

Although  further  theoretical  and  experimental  work  remains  to  be  done  to¬ 
wards  the  analysis  of  accumulated  round-off  error  and  truncation  error  in  the 
new  methods,  it  is  expected  that  they  can  eventually  be  incorporated  into 
efficient  algorithms  for  solving  the  general  equations  of  motion  of  an  earth 


satellite. 
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1.  Introduction 

The  research  reported  in  this  paper  was  motivated  by  a  desire  to  construct 
new  predictor-corrector  methods  for  numerically  solving  the  equations  of  motion 
of  an  earth  satellite. 

Among  many  others,  the  Gauss-Jackson  method  of  solving  these  equations  has 
been  in  use  for  many  years  with  comparatively  good  results.  This  method  is 
based  in  part  on  a  difference  equation  corresponding  to  a  special  differential 
equation  of  the  second  order  (with  first  derivative  absent).  It  appears  to 
owe  to  this  characteristic  some  of  its  advantages  over  other  more  general 
numerical  integration  schemes.  It  is,  moreover,  in  part  a  result  of  applying 
a  particular  process  of  summation  to  the  above-mentioned  difference  equation. 

A  clear  advantage  can  be  shown  [ 5  ]  in  certain  cases  for  the  summed  form  in 
the  control  of  the  propagation  of  round-off  error.  If  arithmetic  is  carried 
out  with  "sufficient"  precision,  the  advantage  is  not  clear.  Neither  is  the 
advantage  of  the  summed  form  clear  if  the  so-called  process  of  double  precision 
accumulation  is  used  with  the  orthodox  form.  This  problem  has  been  very 
briefly  considered  here.  The  results  of  an  experiment  are  given  where  the 
summed  and  unsummed  forms  of  a  simple  difference  equation  are  used  with  and 
without  the  relatively  inexpensive  process  of  double  precision  accumulation. 

With  this  important  exception,  difference  equations  used  in  the  present  are 
unsummed.  We  are  primarily  concerned  with  reducing  local  truncation  (discretiza¬ 
tion)  error. 

The  order  of  accuracy  (Definitions  3.6  and  3.10)  of  the  correcting  difference 
equations  of  the  Gauss-Jackson  method  (summed  or  unsummed) ,  if  kth  differences 
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are  considered,  is  k+1.  In  the  last  decade  much  work  has  been  accomplished 
in  the  analysis  of  linear  multistep  methods  for  solving  differential  equations. 
In  1956,0.  Dahlquist  showed  that  the  order  of  a  stable  k-step  method  could 
not  exceed  k+2  and  that  the  order  could  be  k+2  only  in  a  highly  restricted 
number  of  cases.  Henrici  proved  [5]  this  result  for  a  stable  second  order 
difference  equation,  i.e.,  for  a  difference  equation  corresponding  to  the 
special  equation,  y"=f(x,y). 

If  stability  is  ignored,  a  k-step  linear  difference  equation  corresponding 
to  a  first  order  differential  equation  can  always  be  constructed,  based  on 
the  Hermite  interpolating  polynomial,  with  an  order  of  accuracy  of  2k+l. 

(This  is,  of  course,  as  much  as  one  can  normally  hope  for  in  an  operator  derived 
from  an  oscuiatory  interpolating  polynomial,  because  we  impose  conditions  on 
the  polynomial  and  its  derivative  at  k+1  points.)  Thus  it  is  seen  that  a 
considerable  gap  exists  between  the  maximum  possible  order  of  accuracy  of  a 
k-step  difference  equation  and  the  maximum  possible  order  of  accuracy  of  a 
stable  k-step  difference  equation.  Stability  is  a  necessary  condition  for 
convergence  [5,  p.  217]  of  a  linear  multistap  method.  To  close  the  gap  it  is 
necessary  to  abandon  the  traditional  equally-spaced  multistep  methods. 

In  1964,  Gragg  and  Stetter  introduced  [4]  a  class  of  difference  operators 
for  a  first  order  equation  which  they  called  generalized  multistep  predictor- 
corrector  methods.  The  derivative  value  is  used  at  exactly  one  "non-step" 
point.  These  methods  represent  a  compromise  between  Runga-Kutta  methods  and 
conventional  multistep  methods.  They  use  information  previous  to  the  last 
point  computed,  as  do  t.ne  latter;  and  derivative  evaluations  are  made  part 
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way  through  a  step,  as  in  the  former.  The  authors  prove  that  convergent 

generalized  methods  (of  solving  a  first  order  equation)  with  an  order  of 

accuracy  of  2k+l  exist  for  k  <_  4.  Actually  they  are  able  to  choose  the  non-step 

point  so  that  an  order  of  2k+2  is  obtainable.  Moreover  they  show  that  the 

2k+2 

methods  converge  like  h  uniformly  in  a  given  interval  of  integration.  They 
provide  s  ;ant  computational  evidence  of  the  value  of  their  schemes. 

Butcher  [1],  still  considering  first  order  differential  equations  only, 
developed  an  elegant  method  of  constructing  stable  methods  with  order  of 
accuracy  2k+l  for  k  <  8.  His  method  of  construction  employed  the  familiar 
Hermite  interpolating  polynomial. 

\  In  view  of  the  foregoing  remarks,  it  was  felt  wise  to  studv  multistep 

operators,  generalized  in  the  direction  noted  above,  with  respect  to  a  second 
order  equation  y"=f(x,y).  This  problem  is  interesting  in  itself;  and  it  was 
felt  that  if  satisfactory  results  were  achieved,  a  method  of  this  sort  could 
be  incorporated,  in  analogy  with  the  structure  of  the  Gauss-Jackson  algorithm, 
into  one  which  applied  to  perturbed  satellite  motion  and  could  retain  its 
effectiveness  for  certain  types  of,  if  not  most,  trajectories. 

Our  difference  equations  are  based  on  an  interpolating  polynomial  P(x) 
of  degree  m  such  that 

P(x)  -  g (x)  xeo 

P"(x)  «  g"(x)  XET 

where  o  and  t  are  subsets  of  S^,  a  set  of  k+1  equally  spaced  points.  The  set 
'  oUt  contains  m+1  points,  and  g  is  a  twice  differentiable  function. 
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The  uniqueneaa  of  such  a  polynomial  is  in  question.  It  does  not  depend  on  the 
values  of  g  and  its  second  derivative.  This  problem  of  the  uniqueness  of  the 
generalized  Hermite  interpolating  polynomial,  even  in  its  simplest  form  where 
0“T“Sj^ ,  does  not  appear  to  be  considered  in  the  literature.  Although  we  are 
practically  concerned  only  with  fitting  a  function  and  its  second  derivative  at 
equally  spaced  points,  there  are  many  associated  problems  of  theoretical  interest. 
One  may  replace  the  second  derivative  in  the  above  problem  by  the  n£h  derivative, 
impoae  conditions  on  derivatives  of  mixed  orders,  or  remove  the  requirement 
of  the  equal  spacing  of  the  points. 

A  theory  of  generalized  Hermitian  interpolation  was  needed.  Some  results 
pertaining  to  this  study  are  given  in  Section  2.  Professor  T.  S.  Motzkin  of 
the  University  of  California!  Los  Angeles,  was  consulted  during  the  preparation 
of  this  section.  A  subsequent  paper  by  Professor  Motzkin  and  the  author  will 
treat  this  subject  more  fully. 

The  study  of  conditions  for  the  uniqueness  of  the  generalized  Hermitian 
interpolational  polynomial  leads  to  the  consideration  of  the  question  of 
independence  of  certain  sets  of  polynomials,  to  the  definition  of  generalized 
Vandermonde  determinants  and  to  the  study  of  Hankel  determinants  associated 
with  the  sequence  of  coefficients  in  the  Maclaurin  series  for  logn(l+x) . 

In  Section  3  some  of  the  above  theory  is  used  to  construct  stable  Uigh- 
order  k-step  methods  of  numerical  integration  generalized  as  above.  We  are 
concerned  with  relatively  small  k,  k  <  8.  The  difference  operators  constructed 
depend  on  a  parameter  0.  Graphs  are  included  indicating  the  ranges  of  0  where 
stability  occurs,  in  some  successful  cases.  The  constructive  process  was  used 
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also  to  produce  stable  generalized  explicit  methods  of  order  of  accuracy  2k, 
not  known  to  exist,  for  a  first  order  differential  equation. 

A  final  section  reports  the  results  of  simple  computational  experiments 
comparing  machine  time  required  for  the  new  operators  with  that  required  for 
a  standard  operator  to  compute  a  solution  with  comparable  accuracy.  The 
Appendix  contains  the  coefficients  in  constructed  equations. 

Theorem  2.17  is  due  to  T.  S.  Motzkin.  The  truth  of  Lemma  2.3  was  de¬ 
monstrated  by  R.  B.  Barrar. 
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2.  Generalized  Hermite  Interpolation , 

Our  method  of  constructing  a  high-order  generalized  k-step  difference 
equation  to  solve  numerically  the  equation 

y"  -  f(*,y) 

requires  a  generalization  o'  the  notion  of  Hermitian  interpolation.  The  usual 
k  +  1  point  second-order  Hermite  intarpolational  polynomial  P(x)  is  of  degree 
3k  +  2  and  for  a  twice  differentiable  function  g(x)  satisfies 


(Xj)  =  (xt) 


y  =  0,1,2  i  =  0,1,. ..,k 


g 


In  general,  if  higher  derivatives  appear,  all  derivatives  of  lower  order  are 
considered.  See  fll].  Here  we  wish  to  fit  the  function  and  second  derivative 
and  not  to  impose  conditions  on  the  first  derivative.  We  consider  here  only 
the  case  of  equally  spaced  points  x^.  The  symbol  X  will  always  denote  the 
setfx^Xj,. ..  ,.<k)  . 

The  generalization  of  the  usual  osculatory  interpolation  consists  in 
defining  a  polynomial  P  of  degree  m,  if  possible,  by  imposing  m  +  1  conditions 
of  the  type 

P(P)  (\)  =  g(p)  (xt) 
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with  0  <  p  <  n  and  0  <  i  <  k.  For  any  such  set  of  conditions,  a  unique  poly¬ 
nomial  is  or  is  not  defined.  Conditions  for  the  uniqueness  are  of  interest. 
One  may  begin  with  the  following  determinantal  form  [9], 


p(x) 

1 

X 

2 

X 

•  •  • 

xj 

xj+1 

m 

X 

*<V 

1 

xo 

2 

xo 

•  •  • 

xj 

0 

xj+1 

xo 

m 

xo 

gO^) 

1 

X1 

2 

X1 

•  .  . 

j+1 

x^  ... 

m 

..  xx 

e(xk> 

1 

xk 

2 

*k 

.  . . 

< 

j+1 

*k 

m 

g'  (xQ) 

0 

1 

2xq 

. . . 

.  j-i 
Jxo 

(  j+1)  Xq 

m- 1 

..  mxQ 

g'(xk> 

0 

1 

•  j-i 
Jxk 

(J+l)xjj 

m-1 

..  mx^ 

g(n>  (*0> 

0 

0 

o 

•  •  • 

*  7 

J  • 

(j+1)  IXq 

ml  m-n 

(m-n)!  A0 

g<n>  (xk) 

0 

0 

0 

•  •  • 

ji 

(j+1)  lxk  ... 

m'.  ..m-n 

(m-n) ! 

defining  the  usual  Hermite  interpolating  polynomial  P(x).  Here  m= (n+1) (k+l)-l . 
For  m  <  (n+1) (k+l)-l,  the  more  general  determinantal  form  is  obtained  from  the 
above  form  by  eliminating  rows  corresponding  to  conditions  not  imposed  and 
eliminating  an  equal  number  of  columns  on  the  right.  The  coefficients  of  the 
interpolating  polynomial  P  are  uniquely  determined  if  and  only  if  the  fundamental 


28  October  1966 


10 


SP-2631 


4 


determinant  which  results  from  deleting  the  first  row  and  first  column  from 
the  general  form  (2.1)  is  different  from  zero.  We  call  this  determinant  a  genera¬ 
lized  Vandermonde  determinant.  Some  special  cases  (Hermite  interpolation)  have 
been  carefully  considered.  They  satisfy  the  following  condition: 

If  a  derivative  condition  is  given  at  a  point,  then  conditions 
cn  all  lower  derivatives  are  given  at  that  point. 

It  is  well  known  that  uniqueness  occurs  in  these  cases. 

In  this  paper  we  are  interested  primarily  in  another  special  case,  that 
resulting  from  the  imposition  of  a  condition  on  the  polynomial  and  its  nth 
derivative  at  all  k  4  1  points  (n  <_  k  4  1).  The  resulting  generalized 
Vandermonde  determinant  is 


(2.2)  D(k,n,X)  =  |M(k,n,X) 


2 

2k4l 

1 

xo 

xo 

X0 

2 

2k4l 

1 

X1 

X  ^  •  •  • 

...  x  ^ 

2 

2k4l 

1 

x,  .  .  . 

k 

...  X. 

k 

n 

1-n 

„  2-n 

„  2k4l-n 

Vo 

Vo 

Vo 

2k4l  0 

n 

l-n 

2-n 

2k4l-n 

Vi 

Vi 

a2Xl 

a2k4lXl 

n 

l-n 

c 

1 

2k4l-n 

U  ft. 

a.x. 

1  K 

a.x.  ,  = 

2  k 

<T  v 

~2k4l  k 

I 


) 
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where  =  i(i  -  1)  (i  -  2)...(i  -  n  +  1).  If  M(k,n,X)  is  non-singular,  a 
uniquely  defined  interpolating  polynomial  may  be  associated  with  it.  Denote 
this  polynomial  by  P(k,n,x).  The  following  lemma  is  true  more  generally.  For 
this  report  only  the  case  n  =  2  is  of  interest. 


Lemma  2.  3.  Let  m  =  2k  +  1.  M(k, 2,X)  is  non-singular  if  and  only  if  the  poly¬ 
nomials  (x  4  x0)ra,  (x  4-  xp1", . . . ,  (x  4-  x^)®,  m(m  -  1)  (x  +  x^)1”"2,  m(m  -  1)  (x  4-  x^1" 
...,m(m  -  1)  (x  4-  x^)  are  linearly  independent. 


Proof.  Denote  the  polynomials  as  ordered  by  P^(x)  ,?2(x)  ,P^(x) , . . .  ,Pro+^(x) . 
Expand  P^(x)  in  powers  of  x: 


Pi(x) 


1  =3i4-1 

l 

j=l 


aijX 


,m- j4-i 


i  »  l,  2,. . .  ,m  4-  1 


and  consider  the  matrix  (a^).  For  \  >  1  entries  fll  x+1,  a2  x+1,  • .  •  in  the 
(A4-l)th  column  are: 


*0’  0 


m 


-  «  (V2  5  Xk’2 


Division  now  " 


by  ^ )  gives  the  (\  4-  1)  th  column  of  M(k,2,X).  For  the  first 


two  columns  the  result  is  obvious. 
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Define  the  translation  operator 
polynomials 


on  the  vector  space  of  nth  degree 


(2.4)  T  P(x)  *  P(x+a)  PeV 
a  n 


Obviously  this  linear  translation  is  non-singular.  (Apply  it  to  the  natural 
basis.)  Define  for  real  a: 

X+a  =  {xg+a,  x^+a,  ...»  x^+a} 


Corollary.  jM(k,  2,  X)  j  +  0  <■=■=>  j  M (k r  2,  X+a)  j  i  0. 

Proof.  Because  T  is  non-singular,  the  set  P„ ,  P„,  P  of  the  lemma 

— — —  a  1  i  m+i 

is  linearly  independent  if  and  only  if  the  set  T  P.,  T  P,,  ...,  T  P  1  is 

a  l  a  L  a  m+i 

linearly  independent. 

The  non-singularity  (singularity)  of  M(k,  2,  X)  is  invariant  under  transla¬ 
tion  by  a  real  number.  Because  of  this  property,  we  may  sometimes  write  M(k,  n) 
instead  of  M(k,  n,  X).  It  is  understood  that  there  are  k+1  equally  spaced 
points  and  that  the  rows  of  M(k,n)  are  ordered  to  correspond  to  conditions  on 
the  interpolating  polynomial  and  its  nth  derivative  as  in  (2.2).  The  following 
fact  is  obvious: 

Lemma  2.5.  M(k,  k+1)  Is  non-singular.  M(k,  n)  is  singular  if  n  >  k+1. 

Corollary.  There  exists  a  unique  polynomial  of  degree  2k+l  which  vanishes 
along  with  its  (k+1) th  derivative  at  k+1  equally  spaced  points  for  all  k  >_  0- 
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(2.6) 


can,  by 


(2.7) 


Theorem 


Proof : 
x0’  xl* 


clear 

that  a 

matrix  A 

-  <v 

of  the 

form 

all 

0 

a13 

0 

al5 

0 

0 

a22 

0 

a24 

0 

a26 

a31 

0 

a33 

0 

a35 

0 

0 

i 

a42 

0 

a44 

0 

a40 

elementary  row  operations,  be  triangularized  in  the  usual  manner  giving 


11 


0 

0 


0 

b 

0 

0 


'13 


22 


33 


'15 


24 


35 


44 


26 


46 


2.8.  M(k,  a)  is  singular  if  both  the  following  conditions  are  satisfied 


1)  n  even,  n  >  0 

2)  k  even,  k  0 

Consider  first  the  case  k  =  n  =  2.  We  have  equally  spaced  points 
x^.  By  the  corollary  to  Lemma  2.3  we  may  take  x^=0.  We  use  the 


notation  M  ^  N  to  state  the  fact  that  matrices  M  and  N  are  row-equivalent. 
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Now  it  is 

clear  that 

r 

i 

0 

0 

0 

0 

0 

i 

1 

1 

1 

1 

1 

i 

-1 

1 

-1 

1 

-1 

M(2, 2 ,X)  % 

0 

Q 

2 

0 

0 

0 

0 

0 

2 

6 

12 

20 

0 

0 

2 

-6 

12 

-20 

v_ 

and  that 

r 

i 

0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

1 

i 

0 

1 

0 

1 

0 

(2.9) 

M(2 ,2 ,X)  * 

0 

0 

0 

6 

0 

20 

0 

0 

2 

0 

12 

0 

0 

0 

2 

0 

0 

0 

> 

Triangularizing  the  matrix  found  by  deleting  the  last  row  and  last  column  from 
the  matrix  on  the  right  as  in  the  statement  preceding  the  theorem,  we  have 


r 


M(2,2,X)  ^ 


bn  0 


0 

0 

0 

0 


0 

0 

0 


b13  ° 


b22  ° 


33 


0 

0 

2 


b24  0 


44 


35 


55 


b15  ° 


26 


0 

b 

0 

0 


46 


J 


) 
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No'-'  it  is  apparent  that  M  is  row -equivalent  to  a  triangular  matrix  with  a  zero  on 
the  diagonal,  i.e.,  the  last  row  can  be  annihilated  by  row  operations,  using 
interchanges  if  necessary.  Clearly,  for  even  k,  we  may  take  the  interpolating 
set  X  to  include  the  origin  and  the  remaining  members  of  X  to  be  symmetrically 
distributed  about  the  origin.  If  in  addition  n  is  even,  a  form  analogous  to 
(2.6) can  be  achieved,  i.e.,  for  elements  a^  not  in  the  last  row,  i+j  odd  implies 
a^=0.  The  entries  in  the  final  row  can  then  be  systematically  reduced  to  zero. 

The  above  theorem  shows  that  a  polynomial  of  degree  2k+l  is  not  uniquely 
determined  by  its  values  and  the  values  of  its  nth  derivative  at  k+1  equally 
spaced  points  if  k  is  even  and  n  is  even. 

M(l,2)  is  non-singular  by  Lemma  2.5.  It  can  be  verified  by  hand  computa¬ 
tion  that  M(3,2)  is  also  non-singular.  Machine  computations  of  j M ( k , 2 ) j  for 
k=5,7,9  indicate  that  the  matrices  are  non-singular  in  support  of  the  following 
conj  ecture : 


Conjecture  2.10.  There  exists  a  unique  polynomial  of  degree  2k+l  which  vanishes 
along  with  its  second  derivative  at  k+1  equally  spaced  points  if  k  is  odd. 


A  further  condition  for  the  uniqueness. 

It  has  been  noted  that  a  generalized  Hermitian  interpolating  polynomial 
is  unique  if  and  only  the  associated  generalized  Vandermonde  determinant 
is  different  from  zero.  We  have  also  given  a  condition  for  the  non-singularity 
of  M(k,2)  involving  the  polynomials  {T^x^^+^,  T^x  ^  ^ ;  i=0, 1, 2, . .  .  ,k)  .  We  now 
derive  a  third  necessary  and  sufficient  condition  for  the  uniqueness  of  P(k,2,x). 
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The  linear  transformation  T-T^  on  the  space  of  polynomials  of  degree 
p  has  been  previously  defined  (2.4).  Another  linear  transformation  D=D^,  the 
differentiation  operator,  can  also  be  defined  on  this  space.  See  [10].  Let  T=T 
D  is  uniquely  given  as  a  polynomial  in  T-I. 

D(T-I)  -  (T-I)  -  1/2  (T-I) 2  +  1/3(T-I)3  -  ...  +  -  -  (T-I)p  . 

P 

Lemma  2,11.  If  g(x)  eV^  and  if  g(T)xp»0,  then  g(x)  is  the  zero  polynomial. 

Proof.  The  independence  of  the  polynomials  xP,  TxP,  T2xP,  ...,  TPxP  is  a 
consequence  of  the  fact  that  the  usual  Vandermonde  determinant  is  different 
from  zero.  The  lemma  follows. 

Corollary.  Let  g^(x)eV^  1*0, 1, 2 , . . .  ,p .  Denote  by  S  this  set  of  polynomials. 
Denote  by  S  the  set  g  (T)xp  i*0,2, . . . ,p.  Then  S  is  independent  if  and  only 
if  Sx  is  independent. 

It  is  clear  that  the  nth  order  differentiation  operator  D^  =  Dn  =  Dn(T-I) 
is  given  by  polynomial  multiplication. 

Dn(T-I)  =  d ^ (T-I) 1  . 

i=0  i 

(2) 

For  future  reference  we  note  that  the  sequence  {d;  ;  i=0,l,2,...}  begins 

0  o  i  _i  11  .1  137  _7_ 

u,  u,  x,  ±,  12,  6  ,  180,  10,  ...  . 

We  let  m*2k+l  and  consider  the  question  of  the  independence  of  the 


polynomials 
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D2(T-I)xm 
T  D2(T-I)xm 
T2D2(T-I)xm 

TkD2(T-I)xm 

One  may  write  the  TJ  as  polynomials  in  T-I  and  equivalently  consider  the 
independence  of  the  set 


m 

X 

D2(T-I)xm 

(T-I)xm 

(T-I)D2(T-I)xm 

(T-I)2xm 

(T-I)2D2(T-I)xm 

• 

(T-l)kxm 

• 

(T-I)kD2(T-I)xm 

Finally  using  the  fact  that  (T-I)J=0  for  j  >  m, and  using  the  corollary  to 

Lemma  2.11,  we  consider  the  independence  of  the  polynomials  in  the  operator 

T-I  in  the  latter  set  above.  Writing  A=t-I  we  write  the  matrix  of  coefficient 

2  2  k+1 

of  these  polynomials  relative  to  the  basis  1,  A,  A  t  A  .  Because  the 

above  discussion  is  not  restricted  to  the  second  derivative,  let  n  ^  k+1 
and  write 


m 


T  x 


m 


2  m 
T  x 


Tk  m 
T  x 


(2.12) 


Dn(T-l)  =  d(n)(T-I)n  +  d(^(T-I)n+1  + 
n  n+1 


,(n>  _  .  2k+l 

+  a2k+1u-l) 


The  following  matrix  is  then  written  for  such  n,  (here  n  <  k)  and  is  writte 
(n> 

for  d.  '  for  convenience.  The  composition  of  the  matrix  is  clear.  The 
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k+1  x  k+1  matrix  in  the  lower  left  corner  is  always  the  identity.  That  in 
the  lower  right  is  zero.  The  first  non-zero  element  in  the  first  row  appears 
under  Xn+^,  in  the  second  row  under  Xn+^  \  and  in  the  (k+l)tli  row  under  An  . 
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Denote  this  array  by  M(k,n).  From  the  form  of  the  matrix  of  dimension 
2k+l  x  2k+l  in  the  lower  left  hand  quarter  of  H(k,n),  it  is  evident  chat  the 
latter  matrix  is  non-singular  if  and  only  if  the  matrix  N(k,n), given  by  the 
array  below,  is  non-singular. 


k+1 

A 

k+2 

A 

n+k-2 

•  •  •  A 

n+k-1 

A 

An+k 

^n+k+1 

A2k 

•  •  •  A 

x2k+l 

XkDn(X) 

0 

0 

...  0 

0 

d 

n 

dn+l 

...  dk 

dk+l 

Xk_1Dn(X) 

0 

0 

...  0 

d 

n 

dn+l 

dn+2 

dk+l 

dk+2 

,k-2„n  , 

A  u 

0 

d 

n 

...  d 

n 

dn+l 

dn+2 

dn+3 

dk+2 

dk+3 

Xn+2Dn(X) 

0 

d 

n 

...  d, 
k 

dk+i 

dk+2 

dk« 

*'•  d2k-n 

d2k-n+l 

An_1Dn(X) 

d 

n 

dn+l 

'  ’  *  dk+l 

dk+2 

A 

uk+3 

dk+4 

•"  d2k-n+l 

d2k-n+2 

a  Dn(  x; 

dk 

dk+l 

dn+k-3 

d n+k-1 

d n+k-2 

dn+k 

d2k-l 

d2k 

d"(  X) 

dk+i 

dk+2 

“•  dn+k-2 

d n+k-1 

dn+k 

d n+k-1 

...  d2k 

d2k+l 

Theorem_2 .  M(k,n)  is  non-singular  if  and  only  if  N(k,n)  is  non-singular. 

Some  special  cases  are  given  below  concerning  the  case  of  current  interest, 
i.e.,  n=2.  We  take  all  entries  to  be  positive  and  call  the  resulting  matrix 
N(k,n) .  Clearly  |N(k,2)|  «  (-l)k+1 |N(k,2) |  . 
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N(l,2)  = 


0  1 
1  1 


N(2 , 2)  I  = 


0  1  1 

1  1  11/12 

i  n/12  :/6 


i'N(3,2)  = 


0 

1 

1 

11/12 

1 

1 

11/12 

5/6 

1 

11/12 

5/6 

137/180 

11/12 

5/6 

137/180 

7/10 

=  (1/240) ' 


The  matrix  N(k,n)  can  be  determined  in  a  slightly  different  manner, 
Define  the  sets 


o  =  {1,  X,  X2,  Xk} 

o(n)  =  { D^(X) ,  XDn(X),  X2Dn(X),  XkDn (X) } 


and  let  S  and  Sv  be  the  vector  spaces  generated  over  the  real  numbers  by 

o  and  o(n'  respectively.  The  sets  o  and  are  both  independent  sets  in 

(  ^  \ 

V  .  For  the  set  Cjo  to  be  independent,  the  formula 
dim(S)  +  dim(S(n))  =  dim(S+S(n))  +  dim(SflS(n))  , 


relating  the  dimensions  cf  two  subspaces  to  the  dimensions  of  their  join  and 

(  11  ' 

intersection,  shows  that  it  is  necessary  and  sufficient  that  SOS  '  = 
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A  necessary  and  sufficient  condition  for  this  is  that 
2  k  n 

(a0+ai^+a2^  +•••+3^^  )D  (1)  be  a  kth  degree  polynomial  for  no  kth  deg tee  polynomial 
2  k 

a0+al^+a2^  +-,,+aj^  •  Conditions  on  the  coefficients  a^  give  the  necessary  and 

sufficient  condition  thct  N(k,n)  be  non-singular. 


Associated  with  any  sequence  {a^}  are  the  so-called  Hankel  determinants 


H1  : 
J 


H1  = 


a . 
l 

ai+l 

•  •  •  3  .  ,  ,  .. 

1+j-l 

ai41 

ai+2 

■  .  ■  a 

i+l 

Vj-i 

a 

“i+j 

ai+2j-2 

If  we  write  logn(l+x)  =  'j  a.x1  to  define  the  sequence  {a.}  then  the 

i=0  1  1 

following  is  evidently  true. 


Theorem.  2.16.  M(k,n)  is  non-singular  and  P(k,n,x)  is  uniquely  determined 


if  and  only  if  H^+1  is  non-zero. 


One  can  prove  a  general  theorem  about  Hankel  determinants  which  has  an 
application  in  the  theory  of  generalized  Hermitian  interpolation. 


Theorem.  2.17.  Let  be  the  Hankel  determinants  associated  with  a  sequence 


a0»  a]_»  a2’  '  '  ‘  ’ 

and  H.l  „^0.  Then 
k-2 


Let  k  be  an  integer  greater  than  2  and  assume  H  ■  0 

K—  J. 


F° 

H°  -  V2  ) 2 

Hk  ", „1  .2  (Hk-l) 

Hk-2 
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Proof.  Denote  by  A^  the  array  associated  with  the  Hankel  determinant  H'f. 

-  J  J 

Let  C1  be  the  g,_c_h  column  of  A?"  and  let  A*  be  the  matrix  formed  by  replacing 
'  J  J  *•  J 

the  £th  column  of  a!'  by  1C'!'  ^  .  Let  hf  be  used  to  denote  I  A^  |  . 

—  J  1  J  l  J  1 1  J 

Th?  columns  ^  of  a||  ^  are  ^^near^y  dependent  by  hypothesis.  Thus  the 
equation 


a=k-l. 

<2-18>  Ji  “i  A-i "  * 


may  be  solved  non-trivially .  Because  Hk_2^0  a  solution  to  (2.18)  is  clearly 


(2.19)  ax  =  -Hk_2 


al  £-lHk-2 


£=2,3 . k-1 


Replacing  the  first  column  of  A^  by  this  combination  of  its  first  k-1  columns 
and  computing  in  terms  of  the  entries  of  the  last  row  of  the  resulting  matrix 
and  their  cofactors,  we  have 

V 

hJ=--hh| 
k  a,  k-1 


where 


i=k-l 


-  I 


i=l 


a.  a, 

l  k-2+i 


i=k-l 


i=l 


ak-2+i  iHk-2 


Consider  the  k+1  x  k+1  matrix 


r 


A 


V. 


al  ”•  ak-3  ak-l 


a2  ak-2  ak 


a  ^  ...  a.  1  a.  . 

j  k-1  k+1 


ak-2  ak-l 


a2k-5  a2k-3 
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It  is  clear  that 


H  =  (-l)k+1|A| 


Let  the  columns  of  A  be  called  1=1,2 ,... ,k-l  and  list  the  columns  of 


r\j 


A  A 
’  k-1 

.  A,.  i  for 

IV""  1 

convenience. 

Columns 

of  A: 

C  r 

3°,  *•*' 

k-3°  * 

k-2C  *  k 

Columns 

o£  'STr 

2^  >  ’ 

4°’ 

k-2C  ’ 

c° 

k-1  k-1 

Columns 

of  a^-T 

1C’  2C  * 

3C*  •••• 

k-3°  ’ 

k-2C*  k- 

It  follows  from  these  sequences  of  columns  and  from  the  equations  (2.18),  (2, 
expressing  the  linear  dependence  of  the  columns  of  A^  ,  that 

%  if  air-  !  1 


|A| 


-  ,uw  !!k2  „1 

(  l)  .1  Hk-1 


H 


k-2 


and  hence  that 


H 


h:  - 


0 

K=2.  _  fHl  -v  2 
1  ,2  Hk-1 


^h;.2) 


The  proof  is  complete. 


19) 
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It  is  clear,  of  course,  that  under  the  obvious  hypotheses  we  have  proved 


H 

(2.20)  ^  -  (H”+h 2 


n+1  2  '“k-1 
Hk-2 


*hich  is  really  no  more  general. 
Let  the  sequence  a^,  a^,  a^* 


.  be  defined  by 


log(i+x)  =  l  a£x£  . 


£=0 


From  the  uniqueness  of  the  usual  Hermite  interpolating  polynomial  one  knows 
that  M(k,  1,  X)  is  non-singular  for  all  k.  If  the  H*  are  associated  with  the 

above  sequence,  it  follows  from  Theorem  2.16  that  for  j =1 ,2,3 . 

If  the  sequence  a^,  a^,  a^,  ...  is  defined  by 


log2(l+x)  =  l  a  x£ 
£=0 


and  the  are  associated  with  this  sequence,  then  we  know  from  Theorem  2.8 

and  Theorem  2.15  that  H^=0  for  j  =3 ,5,7 .  We  wish  knowledge  of  for  even 

njyo  for  j=2,4.  These  determinants  appear  as  | N ( 1 , 2 ) |  and  | N (3 , 2 )  |  on  page  21. 
As  a  corollary  to  the  last  theorem  we  have: 

i  2  1 

Corollary.  With  the  last  definition  of  the  H.,  if  H.^O  for  k  >_  4  then  HjVO 

j  K  k 

for  k=2,4,6,...  and  for  such  k,M(k-l,2)  is  non-singular  and  P(k-l,2,x)  is 
uniquely  defined. 
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3 .  Existence  and  Construction  of  Operators 

Consider  a  differentiable  function  y(x)  and  the  differential  equation 

(3.1)  y*  =  g(x,y). 


Most  k-step  methods  of  solving  (3.1)  are  specializations  of  the  formula 


(3.2) 


a.  y  ,+...+  u  +  a„y  =  h(0  g  ,+...  +  0.g  .  +  0„g  ) 

kyv-k  ,y  O'v  k°v-k  l°v-l  0  v 

1  v-1 


where  y,=*y(x.)  and  g.-g(x,,y,).  One  may  associate  with  (3.2)  the  difference 
J  J  111 

operator 


(3.3)  Lh[y(x)]  =  ctky(x-kh)  +  ak_1y(x-(k-l)h)  +  ...  +  a^yCx-h)  +  aQy(x) 


-  h{6ky'(x-kh)  +  B^y  (x-(k-l)h)  +  ...  +  0Jy'(x-h)  +  0oy'(x)}. 


operates  on  any  differentiable  function.  Assume  that  y(x)  possesses  an 
indefinite  number  of  derivatives.  Each  term  in  the  right  member  of  (3.3)  can 
then  be  expanded  by  Taylor's  formula  to  give: 

(3. A)  Lh[y(x)]  =  CQy(x)  +  C^hy^x)  +  ...  +  C^hqy  ^  (x)  +  ... 

The  coefficients  are  given  by 


C„  =  ot„  +  a.  +  .  .  .  +  a 

0  0  1  A 

C1  =  -[a^  +  aa^  +  ...  +  kak  -  (0^  +  0^  +  ...  +  0k) ] 

(3.5) 


C 

q 


=  (-D4[i 


9r  i  /, 


« t  ~ 


1 


+  - 


(q-1)  ! 


(6  +  2q~10, 

1  l 


.+  kq  1Bk)] 

♦ 


for  q=2, 3, .  . . 
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Definition  3.6.  The  difference  operator  (3.3)  is  of  order  of  accuracy  p  if 
C  =C,=...=C  =0  and  C  +1  ^  0. 

0  1  p  p 

The  order  of  a  difference  operator  is  the  first  crude  measure  of  its 
accuracy.  Among  methods  of  the  same  accuracy  p,  the  coefficient  C  which  of 
course  is  independent  of  the  choice  of  function  y(x) ,  may  be  consiu  -red  as  a 
finer  measure  of  the  accuracy.  The  coefficient  must  be  suitably  normalized 
because  it  can  be  made  as  small  as  one  wishes  by  multiplying  (3.2)  by  a 
suitably  small  constant. 

It  is  clear  from  the  definition  of  C_^  that  the  order  of  an  operator  is  p 
if  and  only  if  L^[y(x)]  =  0  whenever  y  is  a  polynomial  of  degree  not  exceeding 
p,  but  is  non-zero  for  some  polynomial  of  degree  p+1. 

There  is  a  theory  [5]  similar  to  the  above  if  one  begins  with  a  special 
second  order  differential  equation  of  the  form 

(3.7)  y"  =  g(x,y) . 

We  give  here  only  the  difference  equation 

? 

(3.8)  a.  y  +  ...  +  ay  ,  +  any  *  h  (6,  g  ,+...  +  B.g  +  8„g  ), 

k'v-k  l'v-1  (Tv  kBv-k  lBv-l  0&v 

the  difference  operator 


(3.9)  L  2  [  y  (x ) )  =  aky(x-kh)  +  a^^^y  (x-(k-l)h)  +  ...  +  a^ytx-h)  +  aQy(x) 
h 


-  h  {Bky"(x-kh)  +  Bk_1y”(x-(k-l)h)  +  ...  +  B^x-h)  +  B0y"(x)}  , 


and  the  corresponding  definition  of  the  order  of  accuracy,  after  expanding,  as  above. 
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Definition  3.10.  The  difference  operator  L  „  is  said  to  be  of  order  of  accuracy 

h 

p  if  Cn  =  C  =  ...  =C  ,  =0  and  C  ^  0  . 

r  0  1  p+1  p+2 


It  is  again  true  that  the  order  can  determined  from  a  knowledge  of  polynomial 

as  a 

finer  gauge  of  accuracy  among  operators  of  order  p. 


annihilated  by  L  „  and  also  that  C  suitably  normalized  here  can  be  used 

h2  p+2 


Method  of  Construction  . 

Associate  with  the  left  side  of  equation  (3.2)  the  polynomial 


(3.11)  p(U  =  ak+1  +  ak?  +  a^2  +  ...  +  aQCk+1 


Definition  3.12.  A  multistep  method  defined  by  (3.2)  is  said  to  be  stable  if 
and  only  if  all  roots  of  p(£)=0  lie  within  or  on  the  unit  circle  and  a  root  of 
modulus  1  has  multiplicity  at  most  1. 

This  definition  applies  to  (3.8)  also.  Here  the  multiplicity  of  a  root 
of  modulus  1  must  not  exceed  2.  A  definition  of  stability  in  general  is  clear. 

Definition  3.13.  A  multistep  method  is  said  to  be  consistent  if  it  has  order 
at  least  1. 

Stability  and  consistency  together  are  necessary  and  sufficient  conditions 
for  methods  defined  by  (3.2)  and  (3.8)  to  converge.  The  concepts  of  stability 
and  consistency  extend  themselves  to  the  slightly  modified  multistep  methods 
considered  here.  We  are  concerned  now  with  the  construction  of  methods  with 
order  of  accuracy  approaching  2k  and  stabLe  in  the  above  sense. 
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In  1956  G.  Dahlquist  proved  (see  [1],  [41,  [5]  for  appropriate  reference) 
that  the  order  of  a  stable  k-step  method  as  defined  by  (3.12)  cannot  exceed 
k+2  and  can  be  equal  to  k+2  in  only  fairly  unusual  circumstances.  Henrici 
extended  this  result  [5]  to  apply  to  methods  for  second  order  equations.  Thus 
a  method  of  order  greater  than  k+2  is  divergent.  Gragg  and  Stetter  [4]  and 
Butcher  [1)  apparently  independently  showed  that  by  including  a  term  6g 

v—  0 

where  x^_0=x^-0h  in  equation  (3.2)  it  is  possible,  for  an  appropriate  value 
of  0  in  the  interval  0  <  0  <  1  and  for  small  k,  to  determine  the  coefficients 
so  that  the  resulting  method  is  stable  and  has  order  of  accuracy  2k+l  (sometimes 
2k+2) .  This  theory  was  applied  only  to  a  method  for  first  order  equations. 

Gragg  and  Stetter  proved  the  assertion  for  k  _<  4,  and  exhibited  some  implicit 
and  explicit  methods  for  such  k.  They  also  proved  that^ under  reasonable 
hypotheses  on  the  accuracy  of  the  starting  values  and  the  differentiability  of 
the  solution  function,  the  pth  order  methods  produced  true  pth  order 
convergence  uniformly  on  a  fixed  interval  of  integration.  Butcher  used  the 
usual  Hermite  interpolation  theory  to  construct  high-orddr  generalized  implicit 
methods  for  first  order  equations  and  for  k  <  8. 

The  present  author  has  used  as  a  point  of  departure  the  work  described 
above.  To  extend  the  theory  to  apply  to  multistep  methods  of  solving  a  special 
second  order  equation  y"*g(x,y)  in  a  constructive  way,  the  Hermite  interpolational 
theory  needed  to  be  generalized  in  a  direction  apparently  not  before  considered. 
Thus  was  partially  accomplished  in  Section  1.  This  extended  theory  enables  us 
to  define  for  k  <  8  high-order  implicit  and  explicit  stable  generalized  k-step 
methods  of  solving  the  second  order  equation  mentioned  above  as  well  as  explicit 
methods  for  3  <  k  <  8  for  first  order  equations. 
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Consider  again  a  set  of  k+1  equally  spaced  points  xv_^»  xv_k+i*  ’***  x 
on  the  line  and  a  sufficiently  often  differentiable  function  f(x).  Under 
certain  conditions  we  have  already  shown  the  existence  of  a  unique  generalized 
Hermitian  interpolating  polynomial  P(x)  =  P(k,n,x)  such  that 

(3.13)  P(x)  =  f (x )  I 

)  X  £  S 

P(n>-£(n)(x)  I 

It  is  possible  to  generalize  still  further.  Let  m  be  an  integer  such 
that  0  _<  m  <  2 k+1.  Let  a ,  and  t  be  subsets  of  with  the  set  out  containing 

m+1  points.  A  unique  polynomial  P(x)  =  P(o,  t,  k,  n,  x)  of  degree  m  again 
exists  such  that 

P(x)  =  f (x) 

(3.14) 

P'n)(x)  =  f(n)(x) 

if  and  only  if  the  fundamental  determinant  is  different  from  zero.  Although 
we  have  not  extended  theoretical  results  to  this  more  general  setting,  computa¬ 
tion  of  the  determinant  by  machine  shows  that  this  uniqueness  "often"  occurs. 

Recall  the  general  determinantal  form  derived  from  (2.1)  by  eliminating 
rows  corresponding  to  conditions  not  imposed  and  by  eliminating  a  corresponding 
number  of  columns  on  the  right,  reducing  the  degree  of  the  polynomial.  The 
determinant  is  an  mth  degree  polynomial,  Q(x)  .  If  the  fundamental  determinant, 
A,  is  not  zero,  Q(x)=0,  and  we  can  compute  it  in  terms  of  the  entries  in  the 
first  column  and  their  cofactors.  Dividing  by  A,  we  have 


x  e  o 

x  e  t 
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(3.15)  P(x)  5  l  hi(x)f(xi)  +  l  hi(x)f(n)(xi)  . 

X^£0  X  £X 

The  functions  h^,  are  polynomials  of  degree  m  and  depend  on  a,  x,  k,  and  n. 
They  generalize  the  usual  Hermitian  interpolational  functions  of  the  first 
and  second  kind  respectively.  For  n=l  and  a=x=Sk  it  is  well  known  that  they 
have  a  simple  analytical  form.  (This  is  true, of  course, also  when  all  derivatives 
less  than  or  equal  to  the  nth  appear  in  a  form  analogous  to  (3.15)).  In  general, 
however,  no  simple  form  was  found  for  the  functions  h..(o,  x,  k,  n,  x)  and 
h^(a,  x,  k,  n,  x) .  For  this  study  they  were  computed  from  the  determinantal  form. 

Most  familiar  methods  of  numerical  integration  are  obtained  as  a  result  of 
polynomial  approximation.  We  note  without  proof  that  by  proper  choice  of 
n,  a,  x,  x  one  may  obtain  from  (3.15)  all  Adams  type  difference  equations  among 
a  variety  of  other  methods.  For  example, by  specializing 

x  =  x 

v 

n  =  1 

o  =  {x  } 

v-1 

X  —  { X  -  ,  X  .  .  .  ,  X  •. 

v-1  v-2  v-kl 

we  have  the  first  order  explicit  Adams  formula- -the  k-step  Adaras-Bashforth 
equation  which  is  usually  derived  by  integrating  both  sides  of  (3.1)  between 
proper  limits. 

Definition:  The  method  defined  by  (3.2)  or  (3.8)  is  called  explicit  if  6q=0 
and  is  called  implicit  if  Bq^O. 
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This  notation  will  also  be  applied  to  methods  involving  a  "non-step" 
point  about  to  be  described. 

The  formulas 


(3.16) 


a,  v  ,  +a,  ,  v  ,  . .+. 
k' v-k  k-1  v-k+1 


gj  =  g(xj.Yj) 


■+Vv  -  hXVk+ Vi Vk«+- • -+ee Ve+lW 


x  A  =  x  -0h  0  <  0  <  1 

v-6  v 


for  numerically  solving  the  equation  y^=g(x,y)  are  similar  to  the  usual 
formulas  with  the  exception  of  the  term  hn8gg^  g  on  the  right.  We  wish  to 
determine  the  coefficients  a..,  8n. »  80  so  that  (3.16)  is  stable  and  has  order 
of  accuracy  approaching  2k  for  k  <  8.  If  (3.16)  is  rewritten 


(3.17)  hnB6u_e  =  l  Vv-i  *  h"  l  BiVl  1 


<  > 


I 


) 


one  sees  that  if  Oq=  1,  the  remaining  coefficients  can  be  determined  from  the 
nd}  derivatives  of  the  Hermitian  interpolational  functions 
n(o,  t,  k,  n,  x) ,  h(o,  t ,  k,  n,  x) .  Before  considering  the  question  of 
stability,  note  that  by  setting  t=S^-x^  one  eliminates  the  condition 


p<V  -  /nVu) 


The  term  involving  the  derivative  at  x^  is  absent  from  (3.16).  Thus  we  may 
choose  Bq=0  and  an  explicit  method  is  constructed.  The  order  of  accuracy  of 
the  method  is  reduced  by  1. 
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It  remains  to  examine  the  question  of  stability  of  the  method.  Writing 
6  as  independent  variable  instead  of  x,  consider  the  equation 
i=k 

(3.18)  l  a,(6K  -  0  0  <  6  <  1  . 

i=0  1 

Let  f^(0),  •••»  be  the  roots  of  (3.18)  and  let  n=2  in  (3.16). 

Consistency  of  the  method  implies  that  (3.18)  always  has  a  double  root 
Define  the  function 

N(e)  =  max  ' £  (6) [  . 

2<i 

Of  interest  are  the  0  such  that  N(0)  <  1.  Motivated  by  the  concept  of 
strong  stability  [5]  only  such  0  are  admitted.  The  process  as  described 
above  will  fail  if  the  generalized  Hennite  interpolating  polynomial 
P  resulting  from  a  definition  of  subsets  o  and  t  is  not  unique.  As  seen, 
this  always  happens  if  k  is  even  and  o=t*S^.  If  a  unique  polynomial  is  not 
defined,  we  may  redefine  o  and  t  and  try  again,  although  we  wish  to  maximize 
the  number  of  conditions  in  order  to  maximize  the  order  of  the  resulting 
operator . 

The  process  also  fails  if 

N(6)  =  N  (0)  >  1  O<0<1. 

o,T  -  ~  “ 


In  this  case  we  may  again  redefine  t  and  try  again. 
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Figures  1  through  5  show  the  function  N(9)  in  some  successful  cases 
where  stable  high-order  operators  were  found.  For  n=ls2  (applying  to  a  first 
or  special  second  order  equation  respectively)  the  curves  are  labelled  by  the 
step  number  k  and  the  set  t=S, -t  which  is  the  set  of  interpolation  points  at 

iC 

which  the  derivative  condition  is  omitted.  For  example  k=4,  t={x  1  labels 

*  v 

the  function  N ( 0 )  associated  with  the  highest  order  explicit  0-family  of 
operators  of  the  type  considered.  Labels  omitting  t  are  meant  to  imply  t=0  . 

Again  it  is  remarked  that  in  most  cases  considered  we  have  relied  on  an  <ip  proximo  tc 
machine  computation  as  a  criterion  for  the  uniqueness  of  P.  In  all  cases  o-S^. 
Although  usually  of  a  lower  order  of  accuracy,  we  have  included  explicit 

methods  in  the  investigation  because  they  require  only  one  prediction,  y  , 

v— 0 

whereas  implicit  methods  require  two,  y^  and  y^_^.  It  was  felt  that  they  could 
be  of  some  practical  value.  We  are  concerned  primarily  with  solving  a  special 
second  order  differential  equation;  however,  we  have  given  some  explicit  metnods 
for  solving  a  first  order  equation  because  they  have  not  been  included  in  the 
literature . 

The  error  constant  C (0 )  associated  with  the  0-family  of  integration  methods  is 
a  smooth  function  of  0  .  If  N ( 0 )  <  1  in  some  interval  0^  <_  0  <_  0^,  it  is 
reasonable  to  attempt  to  find  a  value  of  0  in  the  above  interval  which  minimizes 
Jc|.  That  problem  has  not  been  considered  here  because  the  resulting  reduction 
in  error  would  usually  be  of  a  lower  order  than  the  one  we  are  principally  investi¬ 
gating.  Results  of  computation,  however,  show  that  the  error  constants  C(9) 
associated  with  the  maximal  order  4-step  explicit  0-family  of  operators  for 
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a  special  second  order  differential  equation  (curve  k“4;  Fig.  4)  can  be  made 
to  vanish  in  the  range  of  stability  by  an  appropriate  choice  of  0.  Thus 
another  complete  order  of  accuracy  can  be  attained  in  this  case  (see  Table  1 

below).  We  have  not  determined  that  particular  value  of  0. 

A  tabulation  of  the  maximal  order  of  accuracy  observed  for  stable  k-step 

operators  constructed  in  the  investigation  is  given  below  for  k  <  8.  Although 
column  two  is  well  known  [1],  it  is  included  here  for  completeness. 


TABLE  1 


Orders 

of  Accuracy 

k 

for  y' 
explicit 

=f(x,y) 

implicit 

for  y" 
explicit 

=g(x,y) 

implicit 

1 

2k 

2k+l 

2k-l 

2k 

2 

2k 

2k+l 

2k- 1 

2k-l 

3 

2k 

2k+l 

2k- 1 

2k 

4 

2k 

2k+l 

2  k- 1* 

2  k- 1 

5 

2k 

2k+l 

2k- 2 

2k 

6 

2k-l 

2k+l 

2k-3 

7 

2k+l 

2k- 2 

*  2k  can  probably  be  attained  here. 
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FIGURE  1.  N(6)  FOR  n=l  AND  INDICATED  k,  t  . 
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FIGURE  2.  N(0)  FOR  n=l  AND  INDICATED  k,T. 
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FIGURE  4.  N(0)  FOR  n=2  AND  INDICATED  k,  t. 


FIGURES.  N(6)  FOR  n=2  AND  INDICATED  k 
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This  method  of  constructing  stable  operators  was  suggested  by  the  work  of 
Butcher  [1]  relative  to  the  case  described  earlier,  which  we  can  obtain  by 
specializing  n=l,  o=t=S^.  Although  the  coefficients  in  this  case  are  calculated 
directly  (without  determinants)  from  well  known  formulas,  it  is  not  clear  how 
the  roots  of  (3.18)  are  determined.  In  any  event,  in  our  work  they  are 
calculated  (approximately)  on  computing  machinery;  and  although  great  care 
has  been  taken  to  insure  that  adequate  significance  has  been  retained,  without 
the  inclusion  of  an  error  analysis  of  our  results, we  have  not  given  a  formal 
proof  of  stability. 

The  theory  of  multistep  methods  for  a  special  equation  of  the  second  order 
parallels  that  for  a  first  order  equation.  It  is  certain  that  in  the  former 
case,  also,  a  theorem  can  be  proved  giving  the  uniform  pth  order  convergence 
in  a  fixed  interval  cf  integration  for  our  generalized  methods  and  can  be 
patterned  after  the  analogous  theorem  for  the  traditional  method  (cf.  con¬ 
vergence  theorem  in  [4]).  The  order  of  accuracy  of  a  method  is  meaningful 
only  asymptotically  as  step  size,  h,  approaches  zero.  Because  of  this,  and 
because  we  have  already  on  two  occasions  of  necessity  replaced  formal  proof 
by  a  computational  argument,  we  elect  to  test  on  available  computing  machinery 
in  a  fairly  realistic  context  some  of  the  methods  constructed  and  to  demonstrate 
convergence  and  order  of  convergence  in  this  way. 
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4.  Computational  Experiments. 

Some  of  the  methods  given  in  the  appendix  for  step  number  k^4  were  coded 
in  FORTRAN,  and  computations  were  performed  on  the  IBM  7094  in  order  to  give  a 
feeling  for  their  behavior.  As  a  standard  of  comparison^ the  computations 
were  also  done  using  what  Henrici  [  5  ]  calls  the  "Cowell"  method  of  numerical 
integration.  (This  is  still  another  use  of  the  term.)  A  single  standard 
trajectory  was  selected  for  computation. 

The  oioit  used  is  characterized  by  the  initial  elements  a=1.5, 
i=45°,  ft=M=(jj=0,  and  initial  values  as  given  in  [2,  p.42]  .  No  perturbations 
have  as  yet  been  considered.  The  period  of  the  orbit  is  between  155  and  156 
minutes . 

The  well-known  Cowell  method  was  selected  as  a  standard  of  comparison 
because  it  is  a  near -maximal  order  multistep  method  and  because  It  is  the  unsummed 
form  of  the  Gauss-Jackson  method.  If  no  round-off  error  is  present,  the  summed 
and  un8uramed  difference  equations  will  yield  identical  results.  The  truncation 
errors  are  the  same. 

Aa  noted  earlier  of  our  generalized  methods,  "prediction’’  and  "correction" 
normally  occur  at  different  time  points.  In  the  context  of  the  generalization 
it  seems  more  apt  to  speak  of  the  process  as  consisting  of  one  or  two  "initial" 
methods  and  a  "terminal"  method.  Our  terminal  method  may  be  explicit  (open). 

Such  combinations  have  been  selected  from  methods  appearing  in  the  appendix  and 
will  be  described.  The  order  of  the  combination  we  have  called  the  minimum  of 
the  orders  of  the  components.  The  step  number  of  the  combination  we  have 
called  the  maximum  of  the  step  numbers  of  the  members.  Although  the  advantage 
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of  the  generalized  methods  over  the  usual  multistep  methods  increases  with 
the  step  number,  we  have  chosen  for  the  experiments  a  relatively  small  step 
number  k=4  for  simplicity.  The  orders  of  the  methods  tested  vary  from  4  to  6. 

An  error  term  cculd  be  derived  analytically  for  the  new  methods;  however, 
this  would  be  cumbersome  anu  perhaps  practically  meaningless.  Round-off  error 
effects  would  not  be  included.  Accumulated  error  is  analytically  complex 
as  are  instabilities  induced  by  a  bad  fitting  of  initial  and  terminal  equations. 
For  these  reasons  the  results  of  experiments  are  presented  to  suggest  the 
efficiency  of  the  new  methc  is. 

In  addition  to  the  orbit  described  above,  we  have  tested  two  methods  by 
numerically  solving  the  equation  y"=  -y,  with  y(0)  =  0,y'(0)=l. 

The  methods  used  will  be  described.  Initial  or  predicting  methods  based 
on  the  theory  of  generalized  Hermitian  interpolations  of  Section  2  are  here 
and  in  the  Appendix  called  quasi-Hermitian  methods.  They  are  not  necessarily 
stable.  The  following  five  combinations  are  identified  as  Programs  A,  B,  C,  D 
and  E. 


Program  A:  The  maximal  order  Stormer  methcd  (order  4)  as  predictor;  the 
Cowell  method  (order  5)  as  corrector. 

Program  B:  The  quasi-Hermitian  method  I  (order  6)  as  the  predicting  method; 
the  Cowell  method  (order  5)  as  corrector. 


Program  C:  The  quasi-Hermitian  method  II  (order  6)  as  the  initial  method;  the 

generalized  multistep  explicit  method  I  (k=3,  order  5)  as  terminal. 

Program  D:  Two  initial  methods  required:  The  quasi-Hermitian  methods  I  and  III 
(both  of  order  6)  as  initial  methods;  the  generalized  multistep 
implicit  method  I  (k=3,  order  6)  as  terminal. 

Program  E.  The  quasi-Hermitian  method  III  (order  6)  as  the  initial  method;  the 
generalized  multistep  explicit  method  II  (order  7)  as  terminal. 
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Table  2  gives  the  error  in  the  solution  to  y"=-y  with  above  initial 
conditions  for  indicated  integration  interval  h  at  three  different  values  of  x. 
Error  is  experimental  value  less  true  value. 


TABLE  2 

Comparisons  of  Error  in  Sin  x 

Error  X  10^ 

Program  A 


h 


.01 

0 

.05 

-46 

.10 

-1971 

.20 

-88389 

h 

Error  X  108 

Program  A 

10 

142 

20 

5374 

Error  X  lo' 

Program  A 

h 

10 

40 

20 

1351 

at  x=3  radians 

Program  E 

-1 

0 

1 

90 

at  x=100  radians 

Program  E 

0 

-8 

at  x=200  radians 


0 

-2 
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Figures  6  through  9  compare  errors  in  one  coordinate  of  the  described  orbit 
trajectory  as  computed  by  the  indicated  programs  at  the  stated  times  for  three 
given  values  of  h.  Instead  of  the  error,  we  plot  -log^Q  | error]  so  that  the 
ordinate  shows  the  number  of  decimal  digits  of  accuracy  in  the  solution. 

The  truncation  error  in  any  of  the  methods  is  strongly  oscillatory  with 
the  period  of  the  orbit  and  grows  with  time.  In  Figure  10  we  attempt  to 
eliminate  the  periodicity  and  to  get  a  notion  of  the  asymptotic  behavior  of 
the  propagated  truncation  error  in  programs  E  and  A.  We  have  integrated  to 
a  point  P  on  the  trajectory  (represented  by  the.  zero  abscissa  on  the  graph). 

As  the  integration  proceeded,  the  error  at  P  was  plotted  after  every 
revolution  up  to  nine  revolutions. 

Figure  11  compares  the  error  in  cos  x  in  an  experiment  using  the  summed 

and  unsummed  form  of  the  difference  equation.  Stormer's  method  (k=2)  [5] 

was  used.  The  solution  was  coded  in  single-precision  FORTRAN  and  was  run  on 

* 

the  IBM  7094.  The  two  methods  have  identical  truncation  errors.  The  difference 
is  due  to  differences  in  round-off  error  and  error  propagation.  Henrici 
outlines  [6]  a  statistical  theory  of  round--off  error  which  shows  the  summed 
form  to  be  clearly  superior.  The  theory  is  verified  by  computation.  Figure  12 


shows  the  results  of  the  same  experiment  as  illustrated  by  Figure  11  except 


t-V,,,*-  - -  ~c  5  f ....  _ _ I-.--: — !'•  r  -  _  i 

*  ■  *  1  i  *  '  ‘‘  ‘  -  -  - - | - — —  — uvvul«UJ-U  WiWU  [  4.  )  p  •  ± 


-Lj  is  useu.  mis 


process  is  very  inexpensive  because  the  evaluation  of  the  derivatives,  normally 


by  far  the  greater  part  of  the  work,  is  accomplished  in  single  precision. 


Figure  12  shows  that  the  orthodox  form  of  the  difference  equation  combined  with 


double  precision  accumulation  is  superior  to  the  summed  form  in  a  simple  experiment. 


error! 
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FIGURE  6.  PRECISION  IN  COMPARED  METHODS,  t  =  3min 
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h 

FIGURE  7.  PRECISION  IN  COMPARED  METHODS,  l=9min 


NUMBER  OF  CYCLES  h=3-2"2 

FIGURE  10.  PRECISION  IN  4-STEP  COWELL  AND  GENERALIZED 


error 


2000 
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SUMMED  FORM  OF  DIFFERENCE  EQUATION. 
ACCUMULATION  USED. 
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5.  Conclusions. 

The  purpose  of  the  research,  of  which  this  paper  reports  one  phase,  is 
to  construct  new  algorithms  for  numerically  solving  the  equations  of  motion 
of  an  earth  satellite.  These  algorithms  are  close  to  Gauss-Jackson  algorithms 
in  structure.  Recent  developments  in  the  theory  of  multistep  methods  have 
been  applied  to  increase  the  accuracy  of  component  methods.  It  is  hoped  that 
the  final  algorithms  will  increase  the  ratio  of  the  order  of  accuracy  to  the 
step  number  (the  number  of  differences  used).  In  this  way  the  starting  time 
for  the  muitistep  process  may  remain  unchanged  while  the  interval  of  integra¬ 
tion  is  increased.  If  stability  characteristics  and  round-off  growth  are 
satisfactory  in  the  new  algorithms,  the  total  time  consumed  in  numerical 
integration  may  be  reduced. 

The  increase  in  accura  :y  of  the  applied  generalized  multistep  methods  for 
small  step  number  seems  sufficiently  great  to  justify  the  effort  of  testing 
higher-order  methods  that  have  already  been  constructed,  of  fitting  these 
component  methods  into  a  final  algorithm,  and  of  testing  the  result  in  a 
completely  realistic  context. 

The  generalized  theory  had  not  before  been  applied  to  the  equation 
y"=f(x,y).  The  pertinent  global  convergence  theorem  (cf.  [4]  Theorem  3.1) 
is  undoubtedly  true  but  is  yet  fn  he  pr<w<?d ,  along  with  some  stability  theorems. 

No  work  as  yet  has  been  done  towards  deriving  a  summed  form  for  the  new  difference 
equations.  This  form  could  be  unwieldy.  The  question,  in  general,  of  the  need 
for  and  effectiveness  of  a  summed  form  in  this  work  is  being  seriously  con¬ 
sidered  in  various  installations.  Figures  11  and  12  report  a  rudimentary 
experiment  relating  to  this  question. 
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Uniqueness  theorems  (of  which  Conjecture  ?.10  is  an  example)  are  of 
continuing  interest  in  the  theory  of  quasi-Hermite  osculatory  polynomial 
interpolation;  also  of  importance  is  the  question  of  defining  tha  resulting 
weighting  functions  h^(x),  h^(x),  etc.  by  means  of  simple  formulas.  For 
trajectories  with  significant  perturbations,  it  may  be  well  to  base  methods  on 
a  polynomial  o*  the  form 

i“k  i“k  _  i-k  _ 

P(x)  =  l  hi(x)y(xv_1)  +  l  hi(x)y,(xv_i)  +  J  h  (x)y"(x  ) 

i-0  i-0  i-0 

(  cf.  the  Obrechkoff  method)  with  perhaps  nongrid  points  also  appearing.  A 
maximum  possible  order  of  accuracy  is  more  than  3k.  One  can  also  investigate 
the  feasibility  of  using  least  squares  polynomials  or  functions  other  than 
polynomials  for  approximation  in  orbit  trajectory  computation. 

The  author  plans  also  to  explore  the  possible  definitions  of  a  reasonable 
asymptotic  average  rate  (t-*°)  of  accumulations  of  truncation  error  associated 


with  a  multistep  method. 
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Appendix.  Coefficients  f-,»r  Generalized  Methods 

In  this  Appendix  are  listed  some  of  the  constructed  stable  generalized 
methods  and  some  Hermitian  predicting  equations.  A  few  combinations  of  these 
have  been  tested  on  computing  machinery  and  the  results  reported  in  Section  4. 

All  methods  can  be  characterized  by  specifying  a ,  t,  k,  n,  0.  And  by 
labelling  as  Hermitian  predicting  equation  or  generalized  multistep  method,  a 
and  T  have  been  defined  in  Section  2.  Figures  1  and  2  concern  first  order  equa¬ 
tions  (n  =  1)  and  show  approximate  ranges  of  0  for  which  stable  explicit  methods 
exist.  The  coefficients  of  the  corresponding  difference  equations  are  fairly 
easily  computed  and  are  not  included  here.  In  all  remaining  computational  work 
we  have  taken  n  =  2.  In  the  Appendix  a  -  if  unspecified.  Concerning  the 
generalized  multistep  methods,  for  each  k,  T  pair  usually  only  one  stable  method 
is  listed  (or  tested)  corresponding  to  a  convenient  value  of  0.  A  9-family  of 
methods  always  exists. 

All  predicting  (initial)  methods  are  based  on  the  generalized  Hermitian 
(called  here  "quasi-Hermitian")  interpolating  polynomial  whose  unique  existence 
for  given  k,  n,  a,  T  was  discussed  in  Section  2. 

The  form  of  the  quasi-Hermitian  predicting  equation  is 


yv-9  =  Vv-ktVlyv-k+l+* 


40.  y  ,+ti  (ft  f  .+6.  .f  ..+6,f  .). 

rv-1  VHk  v-k  ^k- 1  v-k+1  K1  v-r 


The  form  of  the  generalized  multistep  method  is 


y  =  a  y  .  +CL  .y  .  , ,+. . .  4a,  y„  .+h^(&  f  .  +&  .f  .  ,.+...+8.f  ,+Bf  ^+BAf 
KJv-k  K-l  v-k+1  rv-1  Kk  v-k  Kk-1  v-k+1  1  v-1  ^  v-0  ^0  v 
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The  coefficients  are  given  with  fifteen  digits.  The  error  in  the  co¬ 
efficients  has  not  yet  been  completely  determined.  For  k  =  3,  4  probably 
all  digits  are  accurate.  For  k  =  7  probably  no  more  than  eight  digits  should 
be  relied  upon.  For  intermediate  k  an  intermediate  number  of  digits  are  valid. 
For  convenience,  we  specify  a  and  t  in  complementary  form  a  =  s^  -  a  and 

t  =  s.  -  r. 
k 

Coefficients  for  generalized  multistep  (implicit)  methods. 

I.  k=3  T-4»  e=.3 

0^  =  .164864864864865-10  pQ  =  . 825825825825825 -lo'2 

a2  =  -.297297297297297  8  =  .140196218627591 

a3  =  -.351351351351351  ^  =  .766087516187516 

P2  =  .415871754107048 
P3  =  . 209376042709376 -10'1 

II.  k=4  i= { x  }  6= . 2 

v-3 

ax  =  .1564622416043223  -10  pQ  =  -.  3027438680537280 -10-1 

a2  =  -.9464169569758374-10'1  p  =.1688144735844922 

a3  =  -.5045838547345016  &  =  .7851898795361902 

a4  =  . 3460313538886233 -10' 1  P2  =  .4800924707582710 

f33-0 


P4  =  .4235373677581134 
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III.  k=5 


t  =  4>  e  =.15 


=  . 167638725140126-10 


=  . 278400595805990 


P0  = 


-.670912228649990*10 


P  =  .201329045334982 


-.  150337671295784*10 


Pi- 


78400  7930456861 


=  .466002632892664 


=  .825862328579257-10 


P2  = 
P3  = 
^4  = 
^5 


324155404103746 


-.494775767295745 


-.134034339375600 


-.373963322694766  *10 


Coefficients  for  generalized  multistep  (explicit)  methods. 


I.  k=3 


T  =  {x  } 


=  .2-10 


=  -.  1  *10 


=  . 124007936507937 


=  . 770833333333333 


=  . 111111111111111 


595238095238095* 10 


II.  k=4 


T  —  {  X  } 

V 


“l  " 


. 205804969722280-10 


a2 

a3  = 
a4 


963457924410107 


247233242848194 


152641470035498 


161595975532325 


=  . 704431783153855 


872587915934802-10 


=  154458906599228 


=  -.951881093873159  -10 
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III.  k=5  T-{x  ,x  0} 

v*  v™2 

0^  =  .307716392780799-10 

a2  =  -. 334316264723224  *10 

a3  *  . 130022320315476*10 

a.  =  . 120385824155265 
4 

a5  =  154610307885770 


IV.  k=7  T“{x  , x  J 

v  v-2 

*  .300166612993176*10 

a2  *  -.313714901027128*10 

a3  =  .683807311516383 

a.  =  .134723867946272*10 
4 

a5  =  116171903379676*10 

a6  =  .222241054902451 
a?  =  .439148682531614*10" 


I-.35 

e0  =  ° 

3  -  . 199342712868376 
PL  *  .542723814092476 
P2  -  0 

f33  =  -.656160625268288 

3,  =  .171431112289051 
4 

35  =  .894415771194486*10" 

=  .36 

p0  =  ° 

3  *  .205438536632750 
PL  =  .539144517864215 
P2  -  0 

33  =  -. 592537575632003 

34  =  .484847995513130 

35  *  -. 200181366173655 

36  =  -.711155078646860*10 
3?  «  -.193932029482787*10 
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Coefficients  for  quasi-Hermitian 

I.  k«4  u=t={x  }  0=0 

v 

a1  =  160000000000000* 102 

a2  =  .  340000000000000* 102 

a3  =  i60oooooooooooo* io2 

a,  =  100000000C00000-102 

4 

II,  k=4  0=T={x  }  Q=.2 

V 

500286720000000*  10 
a2  =  . 124635136000000* 102 
a3  =  -.611842560000000*  10 
=  -.342220800000000 

III.  k=4  o — t ={ x  }  6=.3 

v 

a,  =  193209930000000* 10 

J. 

a2  =  .647195490000000  *  10 

a3  =  -.344761190000000*10 

a,  -  -. 922437000000000 *10-1 
4 

IV.  k=5  o={x  }  T={x  ,x 

v  v  v  ■ 

0^  =  -.  168023255813954 *102 

a2  =  . 425581395348837 *102 

Of3  =  -.  315116279069767  *102 

a.  =  -.  200996777408638  -10 
4 

a5  =  -.802325581395349 


60  SP-2631 

predicting  equations  . 

(3^  “  .266666666666667*10 
P2  -  . 146666666666667* 102 
P3  =  .266666666666667*10 

p4-° 

Pj  =  .129576960000000*10 
P2  =  .560138240000000*10 
P3  =  .  96988 1600CJ.7000 

=  -.194560000000000*  10"^ 

=  . 862500525000000 
P2  =  .303038047500000*10 
P3  =  .433679775000000 
P4  =  -. 721777500000000*10-2 

.4>  9=0 

Pt  =  . 275581395348837  *10 
P2  =  . 146666666666667 *102 
P3  =  -.286046511627907*10 

e4  -  0 

P5  =  -.662494520579550*10 
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I 


V.  k=5  a={x  }  t={x  ,x  .}  0=.3 

v  v’  v-4 

0^  =  -.210714719098474*10  p^  =*  .881649550039971 

a2  =  .826695465383721*10  P2  =•  .303759825000000*10 

a3  =  -.668751562570494*10  P3  =  -.713862014978198 

a.  =  -.729126387971761*10  p,  »  0 

a5  =  -.175047890984738  35  =  -.483931917058517*10 

VI.  k=5  o={x  }  t={x  ,x  ,}  0= . 5 

v  v  v-4 

0^  =  .106292724609375*10  ^  =  .341125488281250 

a2  =  .102539062500000*10  &2  =  .386718750000000 

a3  =  -.180944824218750*10  P3  =  -.504272460937500 

=  -.  782023005506312*10  P4  “  0 

a5  =  -. 698852539062500*10_1  P5  -  -.385435703822545*10 

VII.  k=5  a={x  }  t={x  ,x  1 }  0=.3 

v  v  v-1 

0^  =  . 157175067554755 *102  p  =  0 

a2  =  .  818655257876087-102  p  =  -.  233352230729348  *102 

Q!3  =  -.  189533965786168  *103  83  =  -.  911021071886413  *102 

=  .  753013271876087-102  P4  =  -. 263728213229348 *102 

a5  =  .  176496060554755 *102  85  =  -.862500525000000 


VIII.  k=6  o=t={x  }  0=.15 

v 

=  -. 267354410875809 *102  p  =  .258741730223803*10 

a,  =  -.  542303743902377 *102  82  =  . 353107545202970 *102 
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a3  =  . 163415437280526 -103 

a.  =  -. 529860381340549 -102 
4 

a5  =  -. 278255325746131  -102 
a,  -  -.638051094039407 

O 

IX.  k=6  0=t={xv)  6=.3 

Q£1  =  100919643422326  -102 

a2  =  -. 227114950680780  -102 
a3  =  .665492820631938 *102 
a4  =  -.  202747607559376  -102 
Ct5  =  -.  120925226942311  -102 
(X  =  -.  378539202714476 

X.  k=6  o=t={x  }  6=.35 

v 

a.  =  -.  07953505973247  *10 
a2  =  -. 168428545345370  -102 
a3  =  . 480321543101358  -102 
a.  =  ~.  140846861658563  -102 

4 

a5  =  -.900799716188439  -10 

a,  ==  -.317081388125608 
6 

XI.  k=7  o={  x  }  t  =  { x  , x 

v  V 

-  .  426248672016901  -102 
a2  =  .  106521856117869  -104 
a3  =  . 131573055131036 -104 


P3  *  . 935523697444449  -102 
=  . 356880342431044 *102 
3  =  .251550985402339  -10 
P6  =  . 854784560392494 -10”2 

=  . 128647106629316 '1C 
e2  =  .  146235774698821  -1G2 
P3  =  . 386584562983006  -102 
P4  =  .  154121622333216  *102 
=  . 122101145605511 *10 
P6  =  . 884987329763802  -10'2 

=  .  100266487549777  -10 
P2  =  .  105330757175910  -102 
P3  =  . 280367261702224  -102 
P4  =  .  114420776177617  -102 
P5  =  .954041971289129 
66  =  .836463763365225  -10"2 

L)  9-. 15 

=  0 

»  -. 120432808602138- 103 
P3  =  -. 131699354325931-104 
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a4  =  -.473523475391012  *10 4  308091180980903  -10  4 

a5  =  .  112448958145522 *10 4  ^  140803040314973  *L04 

a6  =  . 111881088447489 *104  ^  =  -.  155735015276831  *103 

a?  =  . 693603082392709  -102  25S741730223803  '10 

XII.  k=7  o={x  }  x=(x  ,x  ,  }  6=.36 

v  v  v-4 

a,  =  -.631355595928052-10  P,  =  ,959963291682580 

j.  1 

a2  =  133597097252952 *102  P2  «  .980492460305444*10 

a3  =  .532090581314935 *102  P3  =  . 229394560344614 *102 

aA  =  -.  344398239116463 *102  B4  =  0 

QL  =  -.197321764953814  pc  =  -.239793128470220  *10 

(X6  =  .219776712103132  *10  P6  =  -.258511564830275  *10" 

a?  =  -. 964138913490856 -10"1  =  . 711986507800205  *10-2 


XIII.  k=8  o=t={x  }  6=.15 

v 

=  -.561616705786098-102 

a2  =  -. 486198844276056 *103 

a3  =  -.  246075559120176*103 

a.  ^  . 158204669038749 *104 
4 

a5  =  -. 251506615194783-103 

a,  ■=  -.484551151396271-103 
6 

a?  =  -. 560120578494670*102 


^  =  .372288312743760*10 
P2  =  .  987049040638957 *102 
P3  =  -637864549912463 
P4  =  .  123026245967712 *104 
P5  =  .  635237310870639  *103 
e6  =>  .  975616454903224- 102 
=  .345201215578264*10 


-.  540791966999407 


Pg  =  .467686974850554*10 
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